bioPIR: Livestock grazing modulation of the Biodiversity - Ecosystem functioning relationships in mountain grasslands

Instituto Pirenaico de Ecología (IPE), CSIC, Zaragoza, Spain

Reset all

Ejecutar esta línea de código en caso de que se quiera limpiar la sesión actual.

#Not run by default
#Run to clean this session
rm(list = ls())

Resources

Trait-based ecology theory - https://www.cambridge.org/core/books/handbook-of-traitbased-ecology/D79AC6C55CA7D3977AD297ED30A38EF0

Categorical TRY data management - https://planthub.idiv.de/jupyter/html/try/try_Species%20tolerance%20to%20frost.html

Session setup

---

Package loading

Cargamos los paquetes necesarios para ejecutar todo el script. Los paquetes se instalan automáticamente si es necesario.

#Indicate the required packages
required_packages <- c(
  'shiny',        #Build interactive apps
  'dplyr',        #Easy grammar for data manipulation
  'tidyr',        #Helpful functions for data manipulation
  'data.table',   #Handle large datasets
  'vegan',        #Useful community ecology functions
  'rtry',         #Functions to manipulate raw TRY data
  'funspace',     #Build, analyze and plot functional traits
  'FD',           #Compute functional diversity indices and CWMs
  'mFD',          #Compute functional diversity indices
  'lme4',         #Performs LMMs
  'lmerTest',     #Allow p-values for lmer models
  'performance',  #Get R2m and R2c from lmms
  'DHARMa',       #Check models assumptions for lmms
  'effects',      #Display effects for linear models
  'piecewiseSEM', #Performs piecewise SEM
  'GGally',       #Package for correlation matrix visualization
  'factoextra',   #Package for PCA visualization
  'ggplot2',      #Draw beautiful plots with ggplot
  'ggtree',       #Phylogenetic tree extension for ggplot2
  'ggfortify',    #PCA extension for ggplot2
  'leaflet',      #Draw beautiful interactive maps
  'gridExtra'     #Arrange multiple plots
)

#Indicate the installed packages
installed_packages <- installed.packages()

#Install and/or load required packages from CRAN
for (pkg in required_packages) {
  if (!pkg %in% installed_packages) {
    install.packages(pkg)
  }
  library(pkg, character.only = TRUE)
}

#Install/or load required packages from github
if (!require('U.PhyloMaker')) install_github("jinyizju/U.PhyloMaker"); library('U.PhyloMaker') #Builds phylogenetic trees

Data loading

Cargamos los datos que se quieren analizar. Hay que cambiar la ruta a los archivos si es necesario.

#Set paths to files
Plot_locations_path <- 'D:/DATOS IPE 2019/Proyecto DIGIPIR/Análisis/Scripts/Datos_bioPIR_plot_locations_R.txt'
Plot_variables_path <- 'D:/DATOS IPE 2019/Proyecto DIGIPIR/Análisis/Scripts/Datos_bioPIR_plot_variables_R.txt'
Plot_sp_matrix_path <- 'D:/DATOS IPE 2019/Proyecto DIGIPIR/Análisis/Scripts/Datos_bioPIR_plot_species_matrix_R.txt'

TRYdata_path <- 'D:/DATOS IPE 2019/Proyecto DIGIPIR/Análisis/TRY/39085.txt'
TRYdata_extra_path <- 'D:/DATOS IPE 2019/Proyecto DIGIPIR/Análisis/TRY/39361.txt'
Trait_matrix_final_path <- 'D:/DATOS IPE 2019/Proyecto DIGIPIR/Análisis/TRY/TRYdata_processed_final.txt'
Trait_matrix_complete_final_path <- 'D:/DATOS IPE 2019/Proyecto DIGIPIR/Análisis/TRY/TRYdata_imputed_final.txt'
Trait_matrix_complete_final_reduced_path <- 'D:/DATOS IPE 2019/Proyecto DIGIPIR/Análisis/TRY/TRYdata_imputed_final_reduced.txt'

#Load data
Plot_locations <- read.table(Plot_locations_path, header=TRUE)
Plot_variables <- read.table(Plot_variables_path, header=TRUE)
Plot_sp_matrix <- read.table(Plot_sp_matrix_path, header=TRUE)

TRYdata <- rtry_import(TRYdata_path)
TRYdata_extra <- rtry_import(TRYdata_extra_path)
Trait_matrix_final <- read.table(Trait_matrix_final_path, header = TRUE)
Trait_matrix_complete_final <- read.table(Trait_matrix_complete_final_path, header = TRUE)
Trait_matrix_complete_final_reduced <- read.table(Trait_matrix_complete_final_reduced_path, header = TRUE)

#Join TRY data
TRYdata_extra$OrigUncertaintyStr <- as.character(TRYdata_extra$OrigUncertaintyStr)
TRYdata<- rtry_bind_row(TRYdata, TRYdata_extra)

Data manipulation

Para facilitar la visualización e interpretación de los datos se reordenan los niveles de pastoreo a Low < MedLow < Med < MedHigh < High.

#Rearrange grazing levels
Plot_locations$Plot<-factor(Plot_locations$Plot, levels=c("Low", "MedLow", "Med", "MedHigh", "High"))
Plot_locations$Plot_rev<-factor(Plot_locations$Plot_rev, levels=c("Low", "MedLow", "Med", "MedHigh", "High"))
Plot_variables$Plot<-factor(Plot_variables$Plot, levels=c("Low", "MedLow", "Med", "MedHigh", "High"))
Plot_variables$Plot_rev<-factor(Plot_variables$Plot_rev, levels=c("Low", "MedLow", "Med", "MedHigh", "High"))

Es necesario homogeneizar la variable continua de presión ganadera. Para ello se calculan z-scores separando por sitio.

#Compute separate z-score for quantitative grazing
Plot_variables <- Plot_variables %>%
  group_by(Location) %>%
  mutate(Grazing_pressure = (Grazing_pressure - mean(Grazing_pressure, na.rm = TRUE)) / 
                      sd(Grazing_pressure, na.rm = TRUE)) %>%
  ungroup()

Además, creamos un nuevo dataset en el que se excluyen los plots Low, MedLow y MedHigh de Festuca paniculata en pirineos.

#Remove Paniculata plots
Plot_variables_s <- Plot_variables %>%
  filter(!(Transect == "Paniculata" & Plot %in% c("Low", "MedLow", "MedHigh")))
Plot_sp_matrix_s <- Plot_sp_matrix[-c(21:32), ]

Exploratory analysis

---

Study area

---

Representamos las localizaciones de los plots de ambas zonas de estudio en un mapa interactivo.

#Define colors for grazing categories
pal <- colorFactor(palette = c("#7bc043A0", "#a1da27A0", "#f0f123A0", "#f5a527A0", "#db2f14A0"), domain = Plot_locations$Plot)

#Define colors and labels for the legend
legend_colors <- pal(unique(Plot_locations$Plot))
legend_labels <- unique(Plot_locations$Plot)

location_map <- leaflet() %>%
  #Set the base map
  addProviderTiles(providers$Esri.WorldTopoMap)  %>% 
  setView(lng = -4.24, lat = 39.9, zoom = 6) %>%
  #Add markers
  addCircleMarkers(
    data = Plot_locations,
    lat=Plot_locations$Lat,
    lng=Plot_locations$Long,
    radius = 3,
    color = ~pal(Plot),
    label=Plot_locations$ID_subplot,
    ) %>%
  #Add a visual scale
  addScaleBar(position = c("topright"), options = scaleBarOptions(imperial = FALSE)) %>%
  #Add a legend to the map
  addLegend(
    position = "bottomright",
    colors = legend_colors,
    labels = legend_labels,
    title = "Grazing category"
    )

#Print the map
location_map

Data exploration

---

Composition

Vamos a ver cómo se diferencia la composición de especies de los plots por sito y nivel de pastoreo.

#Perform NMDS
NMDS<-metaMDS(Plot_sp_matrix_s, distance="bray", autotransform = FALSE)

## Run 0 stress 0.1169238 
## Run 1 stress 0.1228992 
## Run 2 stress 0.1213847 
## Run 3 stress 0.1169246 
## ... Procrustes: rmse 0.0002861545  max resid 0.002869321 
## ... Similar to previous best
## Run 4 stress 0.124104 
## Run 5 stress 0.1240986 
## Run 6 stress 0.1220227 
## Run 7 stress 0.1172638 
## ... Procrustes: rmse 0.003529778  max resid 0.04123828 
## Run 8 stress 0.1264493 
## Run 9 stress 0.1199335 
## Run 10 stress 0.119416 
## Run 11 stress 0.1210976 
## Run 12 stress 0.1207251 
## Run 13 stress 0.1192249 
## Run 14 stress 0.1227323 
## Run 15 stress 0.1238014 
## Run 16 stress 0.1200871 
## Run 17 stress 0.123261 
## Run 18 stress 0.1223663 
## Run 19 stress 0.1201749 
## Run 20 stress 0.1174611 
## *** Best solution repeated 1 times

#Print NMDS results
NMDS

## 
## Call:
## metaMDS(comm = Plot_sp_matrix_s, distance = "bray", autotransform = FALSE) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     Plot_sp_matrix_s 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1169238 
## Stress type 1, weak ties
## Best solution was repeated 1 time in 20 tries
## The best solution was from try 0 (metric scaling or null solution)
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'Plot_sp_matrix_s'

#Check stressplot
stressplot(NMDS)

#Extract parameters
plot_df <- scores(NMDS)$sites %>% 
  as.data.frame()

#Add Location and Grazing pressure
plot_df$Location <- Plot_variables_s$Location
plot_df$Plot_rev <- Plot_variables_s$Plot_rev

#Plot parameters
plot_nmds <- ggplot(plot_df, aes(x = NMDS1, y = NMDS2, color = Plot_rev, shape = Location)) +
  geom_point(size = 3, alpha = 0.8) +
  stat_ellipse(linetype = 2, size = 1) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  annotate("label", x = 1.75, y = 2, label = paste("Stress =", round(NMDS$stress, 3)), size = 4, color   = "black", fill = "white") +
  theme_bw() +
  theme(strip.background=element_blank()) + 
  labs(color="Grazing pressure")

#Print
plot_nmds

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos.

Primary production

Miramos la relación entre distintas variables de producción primaria.

#Plot parameters
fig_biomass <- ggplot(aes(x = Biomass, y = NDVI), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,1000) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Primary production (biomass)", y ="Primary production (NDVI)", color="Grazing pressure") +
  facet_grid(~Location)

fig_cover <- ggplot(aes(x = Plant_cover, y = NDVI), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,1000) +
  xlim(0,100) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Primary production (Plant cover)", y ="Primary production (NDVI)", color="Grazing pressure") +
  facet_grid(~Location)

#Print
grid.arrange(fig_biomass, fig_cover, nrow = 2)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos.

Diversity

Le echamos un vistazo a la relación entre las variables de diversidad y las variables de producción primaria.

#Plot parameters
fig_ndvi_rich <- ggplot(aes(x = Richness, y = NDVI), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,1000) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Species diversity (richness)", y ="Primary production (NDVI)", color="Grazing pressure") +
  facet_grid(~Location)

fig_ndvi_shann <- ggplot(aes(x = EShannon, y = NDVI), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,1000) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Species diversity (hill numbers)", y ="Primary production (NDVI)", color="Grazing pressure") +
  facet_grid(~Location)

fig_biomass_rich <- ggplot(aes(x = Richness, y = Biomass), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Species diversity (richness)", y ="Primary production (biomass)", color="Grazing pressure") +
  facet_grid(~Location)

fig_biomass_shann <- ggplot(aes(x = EShannon, y = Biomass), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Species diversity (hill numbers)", y ="Primary production (biomass)", color="Grazing pressure") +
  facet_grid(~Location)

fig_cover_rich <- ggplot(aes(x = Richness, y = Plant_cover), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,100) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Species diversity (richness)", y ="Primary production (plant cover)", color="Grazing pressure") +
  facet_grid(~Location)

fig_cover_shann <- ggplot(aes(x = EShannon, y = Plant_cover), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,100) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Species diversity (hill numbers)", y ="Primary production (plant cover)", color="Grazing pressure") +
  facet_grid(~Location)

#Print
grid.arrange(fig_ndvi_rich, fig_ndvi_shann, nrow = 2)

grid.arrange(fig_biomass_rich, fig_biomass_shann, nrow = 2)

grid.arrange(fig_cover_rich, fig_cover_shann, nrow = 2)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos.

Grazing level

En este apartado se explora la relación entre el nivel de pastoreo y diferentes variables de diversidad, producción primaria y fertilidad del suelo.

#Plot parameters
fig_ndvi <- ggplot(aes(x = Plot_rev, y = NDVI, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,1000) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Primary production (NDVI)") +
  facet_grid(~Location)

fig_biomass <- ggplot(aes(x = Plot_rev, y = Biomass, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,300) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Primary production (biomass)") +
  facet_grid(~Location)

fig_cover <- ggplot(aes(x = Plot_rev, y = Plant_cover, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,100) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Primary production (plant cover)") +
  facet_grid(~Location)

#Print
grid.arrange(fig_ndvi, fig_biomass, fig_cover, nrow = 3)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos.

fig_rich <- ggplot(aes(x = Plot_rev, y = Richness, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Species diversity (richness)") +
  facet_grid(~Location)

fig_eshannon <- ggplot(aes(x = Plot_rev, y = EShannon, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Species diversity (hill numbers)") +
  facet_grid(~Location)

#Print
grid.arrange(fig_rich, fig_eshannon, nrow = 2)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos.

fig_pH <- ggplot(aes(x = Plot_rev, y = pH, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (pH)") +
  facet_grid(~Location)

fig_EC <- ggplot(aes(x = Plot_rev, y = EC, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (EC)") +
  facet_grid(~Location)

fig_OM <- ggplot(aes(x = Plot_rev, y = MO, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (OM)") +
  facet_grid(~Location)

fig_C <- ggplot(aes(x = Plot_rev, y = TC, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (TC)") +
  facet_grid(~Location)

fig_N <- ggplot(aes(x = Plot_rev, y = TN, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (TN)") +
  facet_grid(~Location)

fig_CN <- ggplot(aes(x = Plot_rev, y = CN, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,40) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (CN)") +
  facet_grid(~Location)

fig_P <- ggplot(aes(x = Plot_rev, y = P_Olsen, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (Available P)") +
  facet_grid(~Location)

fig_NP <- ggplot(aes(x = Plot_rev, y = NP, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (ln(NP))") +
  facet_grid(~Location)

#Print
grid.arrange(fig_pH, fig_EC, fig_OM, fig_C, fig_N, fig_CN, fig_P, fig_NP, nrow = 8)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos.

Grazing level (quantitative)

Se vuelve a explorar la relación entre el nivel de pastoreo y diferentes variables de diversidad, producción primaria y fertilidad del suelo. En esta ocasión, se usa una variable cuantitativa de presión ganadera.

#Plot parameters
fig_ndvi_q <- ggplot(aes(x = Grazing_pressure, y = NDVI), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Primary production (NDVI)", color="Grazing pressure") +
  facet_grid(~Location)

fig_biomass_q <- ggplot(aes(x = Grazing_pressure, y = Biomass), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Primary production (biomass)", color="Grazing pressure") +
  facet_grid(~Location)

fig_cover_q <- ggplot(aes(x = Grazing_pressure, y = Plant_cover), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Primary production (plant cover)", color="Grazing pressure") +
  facet_grid(~Location)

#Print
grid.arrange(fig_ndvi_q, fig_biomass_q, fig_cover_q, nrow = 3)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos.

#Plot parameters
fig_rich_q <- ggplot(aes(x = Grazing_pressure, y = Richness), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x="Grazing pressure", y="Species diversity (richness)", color="Grazing pressure") +
  facet_grid(~Location)

fig_eshannon_q <- ggplot(aes(x = Grazing_pressure, y = EShannon), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x="Grazing pressure", y="Species diversity (hill numbers)", color="Grazing pressure") +
  facet_grid(~Location)

#Print
grid.arrange(fig_rich_q, fig_eshannon_q, nrow = 2)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos.

#Plot parameters
fig_pH_q <- ggplot(aes(x = Grazing_pressure, y = pH), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (pH)", color="Grazing pressure") +
  facet_grid(~Location)

fig_EC_q <- ggplot(aes(x = Grazing_pressure, y = EC), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (EC)", color="Grazing pressure") +
  facet_grid(~Location)

fig_OM_q <- ggplot(aes(x = Grazing_pressure, y = MO), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (OM)", color="Grazing pressure") +
  facet_grid(~Location)

fig_C_q <- ggplot(aes(x = Grazing_pressure, y = TC), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (TC)", color="Grazing pressure") +
  facet_grid(~Location)

fig_N_q <- ggplot(aes(x = Grazing_pressure, y = TN), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (TN)", color="Grazing pressure") +
  facet_grid(~Location)

fig_CN_q <- ggplot(aes(x = Grazing_pressure, y = CN), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (CN)", color="Grazing pressure") +
  facet_grid(~Location)

fig_P_q <- ggplot(aes(x = Grazing_pressure, y = P_Olsen), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (Available P)", color="Grazing pressure") +
  facet_grid(~Location)

fig_NP_q <- ggplot(aes(x = Grazing_pressure, y = NP), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Soil fertility (ln(NP))", color="Grazing pressure") +
  facet_grid(~Location)

#Print
grid.arrange(fig_pH_q, fig_EC_q, fig_OM_q, fig_C_q, fig_N_q, fig_CN_q, fig_P_q, fig_NP_q, nrow = 8)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos.

Soil fertility

Le echamos un vistazo a la relación entre variables de fertilidad del suelo y variables de diversidad y producción primaria.

#Plot parameters
fig_ndvi_cn <- ggplot(aes(x = CN, y = NDVI), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,1000) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (CN)", y ="Primary production (NDVI)", color="Grazing pressure") +
  facet_grid(~Location)

fig_biomass_cn <- ggplot(aes(x = CN, y = Biomass), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (CN)", y ="Primary production (biomass)", color="Grazing pressure") +
  facet_grid(~Location)

fig_cover_cn <- ggplot(aes(x = CN, y = Plant_cover), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,100) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (CN)", y ="Primary production (plant cover)", color="Grazing pressure") +
  facet_grid(~Location)

fig_richness_cn <- ggplot(aes(x = CN, y = Richness), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (CN)", y ="Species diversity (richness)", color="Grazing pressure") +
  facet_grid(~Location)

fig_eshannon_cn <- ggplot(aes(x = CN, y = EShannon), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (CN)", y ="Species diversity (hill numbers)", color="Grazing pressure") +
  facet_grid(~Location)

fig_ndvi_p <- ggplot(aes(x = NP, y = NDVI), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,1000) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (ln(NP))", y ="Primary production (NDVI)", color="Grazing pressure") +
  facet_grid(~Location)

fig_biomass_p <- ggplot(aes(x = NP, y = Biomass), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (ln(NP))", y ="Primary production (biomass)", color="Grazing pressure") +
  facet_grid(~Location)

fig_cover_p <- ggplot(aes(x = NP, y = Plant_cover), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,100) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (ln(NP))", y ="Primary production (plant cover)", color="Grazing pressure") +
  facet_grid(~Location)

fig_richness_p <- ggplot(aes(x = NP, y = Richness), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (ln(NP))", y ="Species diversity (richness)", color="Grazing pressure") +
  facet_grid(~Location)

fig_eshannon_p <- ggplot(aes(x = NP, y = EShannon), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (ln(NP))", y ="Species diversity (hill numbers)", color="Grazing pressure") +
  facet_grid(~Location)

#Print
grid.arrange(fig_ndvi_cn, fig_biomass_cn, fig_cover_cn, fig_richness_cn, fig_eshannon_cn, nrow = 5)

grid.arrange(fig_ndvi_p, fig_biomass_p, fig_cover_p, fig_richness_p, fig_eshannon_p, nrow = 5)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos.

Probamos si las propiedades fisico-químicas del suelo se podrían resumir en una única variable.

#Run PCA
pca<-prcomp(na.omit(Plot_variables_s[c(15:23)]), center=TRUE, scale=TRUE)

#Filter df data to match pca rownumber
Plot_variables_s_clean <- Plot_variables_s %>%
  drop_na(15:23)

#PCA plot
autoplot(pca, data = Plot_variables_s_clean, size = 3, colour = 'Plot_rev', shape = 'Location', na.action = na.omit,
         loadings = TRUE, loadings.colour = 'grey25', loadings.label = TRUE, loadings.label.size = 3, 
         loadings.label.colour = 'grey25') +
          theme_bw()

Nota: Parece que con usar el TN o la relación CN sería suficiente.

Functional traits

---

Numerical traits

---

En este apartado procesamos los rasgos funcionales continuos solicitados a TRY para empezar a obtener la matriz final de traits x species.

Data exploration

#Group the input data based on AccSpeciesID, AccSpeciesName, TraitID, TraitName, DataID and DataName and sort by TraitID
#Note: For TraitID == "NA", meaning that entry is an ancillary data
TRYdata_explore <- rtry_explore(TRYdata,
  AccSpeciesID,
  AccSpeciesName,
  TraitID,
  TraitName,
  DataID,
  DataName,
  sortBy = TraitID
)

TRYdata_explore

#Explore the unique values of the OrigValueStr within the selected data
TRYdata_trait_selection <- rtry_select_row(TRYdata, TraitID %in% 146)
TRYdata_values_explore <- rtry_explore(TRYdata_trait_selection,
  AccSpeciesID,
  AccSpeciesName,
  TraitID,
  TraitName,
  DataID,
  DataName,
  OrigValueStr,
  StdValue,
  ErrorRisk,
  Comment,
  sortBy = AccSpeciesName
)

TRYdata_values_explore

Records cleaning

#Remove undesired traits from the trait list
undesired_TraitID <- c(37, 42, 152, 341, 819)
TRYdata_filtered <- rtry_exclude(TRYdata,
  (TraitID %in% undesired_TraitID),
  baseOn = ObsDataID)

#Remove undesired DataID from the trait list
undesired_DataID <- c(1323, 2310, 2385, 7210, 2386, 7211, 2528, 4027, 4028, 7218, 7219, 2311, 19, 20, 3645, 3646, 3647, 4029, 4030, 4142)
TRYdata_filtered <- rtry_exclude(TRYdata_filtered,
  (DataID %in% undesired_DataID),
  baseOn = ObsDataID)

#Exclude outliers
#Identify trait records where the ErrorRisk is larger than 4
#Exclude the records, not the whole observation
#Note: Might result in data losses
TRYdata_exclude <- rtry_exclude(TRYdata_filtered, ErrorRisk > 4, baseOn = ObsDataID)

#Remove duplicates
#Note: Might result in data losses
TRYdata_removed <- rtry_remove_dup(TRYdata_exclude, showOverview = TRUE)

#Explore the imported data
#Select the rows where TraitID is 146, i.e. the data containing the plant leaf CN ratio
#Explore the unique values of the OrigValueStr within the selected data
TRYdata_trait_selection <- rtry_select_row(TRYdata_removed, TraitID %in% 146)
TRYdata_values_explore <- rtry_explore(TRYdata_trait_selection,
  AccSpeciesID,
  AccSpeciesName,
  TraitID,
  TraitName,
  DataID,
  DataName,
  OrigValueStr,
  StdValue,
  ErrorRisk,
  Comment,
  sortBy = AccSpeciesName
)

TRYdata_values_explore

Wide formatting

A partir de aquí se continúa únicamente con los traits con valores númericos continuos. El resto de traits hay que completarlos “a mano” en pasos posteriores.

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYdata_removed, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Extract latitude (DataID 59) and longitude (DataID 60) of each observation within TRYdata
#georef <- rtry_select_anc(TRYdata_removed, 59, 60, showOverview = FALSE)

#Merge the relevant data frames based on the ObservationID using rtry_join_left()
#num_traits_georef <- rtry_join_left(num_traits, georef, baseOn = ObservationID)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

num_traits_wider

Trait matrix

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_t <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}"))) %>%
  select(-"59_Plant lifespan (longevity)_year") # no está bien codificado en los valores continuos

result_t

#Export data
#rtry_export(result_t, file.path('D:/DATOS IPE 2019/Proyecto DIGIPIR/Análisis/TRY', "TRYdata_preprocessed_clean.csv"))

Discrete traits

---

En este apartado continuamos procesando los rasgos funcionales categóricos para ir completando la matriz de traits x species.

N fixation capacity

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Plant nitrogen(N) fixation capacity"]

# extract original data strings and change to lowercase to ease later classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

# define the search strings
searchNames <- c(
    "^no?$|^0$|not an n fixer|non fixer|low|none|not n2 fixing|r",
    "^y(es)?$|^1$|medium|rhizobia|present|frankia|high|^n2-? ?fixing|nostocaceae|^n2?(-| )fixer"
)

# search for the strings defined before
searchResults <- sapply(searchNames, grepl, oriVals)

# name columns of searchResults matrix like corrected searchNames
colnames(searchResults) <- c(0, 1)

# show the number of matches to each category
colSums(searchResults)

# show the original entries for which no match was retrieved
oriVals[rowSums(searchResults) < 1]

# show the number of entries that weren't matched to any category
sum(rowSums(searchResults) < 1)

# show the number of entries that were matched to more that one category
sum(rowSums(searchResults) > 1)

searchResults[rowSums(searchResults) > 1, ] <- FALSE

# get the new values
newVals <- sapply(seq_len(nrow(searchResults)), function(x) {
    paste(colnames(searchResults)[searchResults[x, ]], collapse = ",")
})
newVals[newVals == ""] <- NA

# integrate into TRYSubset
TRYSubset[TraitName == "Plant nitrogen(N) fixation capacity", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_N_fix <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_N_fix

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- result_t %>%
  full_join(result_N_fix, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Photosinthesis pathway

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Photosynthesis pathway"]

# extract original data strings and change all to lowercase to ease later classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

# define the search strings
searchNames <- c("c3", "c4", "cam")

# search for the strings defined before
searchResults <- sapply(searchNames, grepl, oriVals)

# name columns of searchResults matrix like corrected searchNames
colnames(searchResults) <- c(1, 2, 3)

# show the number of matches to each category
colSums(searchResults)

# show the original entries for which no match was retrieved
oriVals[rowSums(searchResults) < 1]

# show the number of entries that weren't matched to any category
sum(rowSums(searchResults) < 1)

# show the number of entries that were matched to more that one category
sum(rowSums(searchResults) > 1)

searchResults[rowSums(searchResults) > 1, ] <- FALSE

# get the new values
newVals <- sapply(seq_len(nrow(searchResults)), function(x) {
    paste(colnames(searchResults)[searchResults[x, ]], collapse = ",")
})
newVals[newVals == ""] <- NA

# integrate into TRYSubset
TRYSubset[TraitName == "Photosynthesis pathway", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_photo_pathway <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_photo_pathway

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_photo_pathway, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Dispersal syndrome

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Dispersal syndrome"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
#valueOverview[order(valueOverview)]

# decode coded entries
strsplit(TRYSubset[DatasetID == 474]$Comment[1], split = ";|,")
codes <- c("G", "W", "H", "B", "M", "N", "P", "O", "Z")
repls <- c("baro", "anemo", "hydro", "ballisto", "myrmeco", "endozoo", "epizoo", "zoo", "zoo")
oriVals[TRYSubset$DatasetID == 474] <- toupper(oriVals[TRYSubset$DatasetID == 474])
for (i in seq_along(codes)) {
    oriVals[TRYSubset$DatasetID == 474] <-
        gsub(codes[i], paste0(" ", repls[i], "chor "), oriVals[TRYSubset$DatasetID == 474])
}
oriVals[TRYSubset$DatasetID == 319 & TRYSubset$OrigValueStr == "yes"] <-
    TRYSubset[DatasetID == 319 & OrigValueStr == "yes"]$OriglName
oriVals[TRYSubset$DatasetID == 351 & TRYSubset$OrigValueStr == "1"] <-
    TRYSubset[DatasetID == 351 & OrigValueStr == "1"]$OriglName
oriVals[TRYSubset$DatasetID == 435 & TRYSubset$OrigValueStr == "1"] <-
    TRYSubset[DatasetID == 435 & OrigValueStr == "1"]$OriglName

oriVals[!grepl("[[:lower:]]", oriVals)] <- NA

# define the search strings
searchNames <- c(
    "herpochor|tumbling|dispersal no|autochor|unassisted|unspecialised|drop to the ground|unspecialised|baro|gravity", # nolint: line_length_linter.
    "^ane$|meteorochor|wind|anemo|boleochor|chamaechor|semachor",
    "eaten|endo-?zoochor|^endo?$",
    "(ecto|epi)-?zoochor|adhesion|^epi$",
    "dysochor|hoarding|squirrel|shrew",
    "(^|[^A-Za-z])ant|myrme(co)?chor",
    "nautochor|water|hydrochor|bythisochor|floating",
    "ballochor|ballistochor|explosive|blastochor|ballistic",
    "([^Aa]|^)biotic|zoochor|grouse|animal|cattle|ox|sheep|bird|deer|horse|pig|rabbit|ma(m)?mal|mouse|goat|hare|donkey|vertebrate|ornithochor|chamois|marten|pacarana|boar|roe|fox|dog|snail|earthworm|badger|trichoptera|monkey|marmot|reptile|fish|civet|beetle|chimpanzee|bat|turtle", # nolint: line_length_linter.
    "(^man)|agochor|human|ethelochor|commerce|vehicle|speirochor|machine|footwear|hemerochor",
    "ombrochor|rain"
)

# search for the strings defined before
searchResults <- sapply(searchNames, grepl, oriVals, ignore.case = TRUE)

# name columns of searchResults matrix like corrected searchNames
# "autochorous", "anemochorous", "endozoochorous", "epizoochorous", "dysochorous", "myrmechorous", "hydrochorous",          "ballochorous", "zoochorous", "anthropochorous", "ombrochorous"
colnames(searchResults) <- c(0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0)

# show the number of matches to each category
colSums(searchResults)

# show the original entries for which no match was retrieved
#oriVals[rowSums(searchResults) < 1]

# show the number of entries that weren't matched to any category
sum(rowSums(searchResults) < 1)

# show the number of entries that were matched to more that one category
sum(rowSums(searchResults) > 1)

# use the searchResults matrix to create new value strings by concatenating all data found
newVals <- sapply(seq_len(nrow(searchResults)), function(x) {
    paste(colnames(searchResults)[searchResults[x, ]], collapse = ",")
})
newVals[newVals == ""] <- NA

# Transformar newVals: si contiene "1", asignar 1; en otros casos, asignar 0
newVals <- ifelse(grepl("1", newVals), 1, newVals)
newVals <- ifelse(grepl("0", newVals), 0, newVals)

# Convertir a numérico
newVals <- as.numeric(newVals)

# integrate into TRYSubset
TRYSubset[TraitName == "Dispersal syndrome", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_dis_syndrome <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_dis_syndrome

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_dis_syndrome, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Drought tolerance

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Species tolerance to drought" & DataID == 34]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

# remove purely numeric values
oriVals[!grepl("[[:lower:]]", oriVals)] <- NA

# define the search strings
searchNames <- c("none|dies",
                 "low|wilts",
                 "medium|recovers",
                 "high|resists")

# search for the strings defined before
searchResults <- sapply(searchNames, grepl, oriVals)

# name columns of searchResults matrix like corrected searchNames
colnames(searchResults) <- c(1, 2, 3, 4)

# show the number of matches to each category
colSums(searchResults)

# show the original entries for which no match was retrieved
oriVals[rowSums(searchResults) < 1]

# show the number of entries that weren't matched to any category
sum(rowSums(searchResults) < 1)

# show the number of entries that were matched to more that one category
sum(rowSums(searchResults) > 1)

# use the searchResults matrix to create new value strings by concatenating all data found
newVals <- sapply(seq_len(nrow(searchResults)), function(x) {
    paste(colnames(searchResults)[searchResults[x, ]], collapse = ",")
})
newVals[newVals == ""] <- NA

# Transformar newVals
newVals <- ifelse(grepl("3", newVals), 3, newVals)
newVals <- ifelse(grepl("2", newVals), 2, newVals)

#-------------------------------------------------------------------------------

# Select data of interest
TRYSubset2 <- TRYdata_removed[TraitName == "Species tolerance to drought" & DataID == 34]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset2$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

# remove non-numeric values
oriVals[!grepl("\\d", oriVals)] <- NA
oriVals <- as.numeric(oriVals)

# reclass values
oriVals <- oriVals + 1 #originals from 0-5 to 1-6
breaks <- seq(1, 6, by = 1.25)
labels <- c("1", "2", "3", "4")
oriVals_class <- cut(oriVals, breaks = breaks, labels = labels, include.lowest = TRUE)

#Join into TRYSubset
newVals = coalesce(as.numeric(newVals), as.numeric(oriVals_class))
TRYSubset[TraitName == "Species tolerance to drought", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_drought <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_drought

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_drought, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Frost tolerance

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Species tolerance to frost"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

# remove purely numeric values
oriVals[!grepl("[[:lower:]]", oriVals)] <- NA

# define the search strings
searchNames <- c("sensitive",
                 "high|resistant|freezingexposed")

# search for the strings defined before
searchResults <- sapply(searchNames, grepl, oriVals)

# name columns of searchResults matrix like corrected searchNames
colnames(searchResults) <- c(0, 1)

# show the number of matches to each category
colSums(searchResults)

# show the original entries for which no match was retrieved
oriVals[rowSums(searchResults) < 1]

# show the number of entries that weren't matched to any category
sum(rowSums(searchResults) < 1)

# show the number of entries that were matched to more that one category
sum(rowSums(searchResults) > 1)

# use the searchResults matrix to create new value strings by concatenating all data found
newVals <- sapply(seq_len(nrow(searchResults)), function(x) {
    paste(colnames(searchResults)[searchResults[x, ]], collapse = ",")
})
newVals[newVals == ""] <- NA

#-------------------------------------------------------------------------------

# Select data of interest
TRYSubset2 <- TRYdata_removed[TraitName == "Species tolerance to frost"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset2$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

# remove non-numeric values
oriVals[!grepl("\\d", oriVals)] <- NA
oriVals <- as.numeric(oriVals)

# convert farenheit into °C
oriVals[TRYSubset2$OrigUnitStr == "farenheit"] <- 5 / 9 * oriVals[TRYSubset2$OrigUnitStr == "farenheit"] - 32

# remove Frost hardiness trait data
oriVals[grepl("Frost Free Days, Minimum", TRYSubset2$DataName)] <- NA

# reclass values
oriVals_class <- ifelse(oriVals < 0, 1, 0)

#Join into TRYSubset
newVals = coalesce(as.numeric(newVals), as.numeric(oriVals_class))
TRYSubset[TraitName == "Species tolerance to frost", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_frost <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_frost

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_frost, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Plant woodiness

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Plant woodiness"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

oriVals[TRYSubset$DatasetID %in% c("9", "86") & oriVals == "0"] <- "non-woody"
oriVals[TRYSubset$DatasetID %in% c("9", "86") & oriVals == "1"] <- "semi-woody"
oriVals[TRYSubset$DatasetID %in% c("9", "86") & oriVals == "2"] <- "woody"
oriVals[TRYSubset$DatasetID %in% c("9", "86") & oriVals == "3"] <- "woody"
oriVals[TRYSubset$DatasetID %in% c("73", "118", "153", "699") & oriVals == "0"] <- "non-woody"
oriVals[TRYSubset$DatasetID %in% c("73", "118", "153", "699") & oriVals == "1"] <- "woody"

# repair hemi-epiphyte matching problem
oriVals <- gsub("emi-epiph", "emiepiph", oriVals)

# remove purely numeric values and others that have no lowercase character included
oriVals[!grepl("[[:lower:]]", oriVals)] <- NA

# define the search strings
searchNames <- c(
    "(nono?-?woody|fibrous|^0$|false|^none?$|no(t|n) woody)",
    "woody? at base|woody rootstock|woody base|semi-woody|basal",
    "(?<!emi)epiph(y|i)tes?",
    "(h|s)emi-?epiph(y|i)tes?",
    "liana",
    "graminoid|grass",
    "^w$|true|^y$|suffrutex|^1$|(^|/)woody(/|$)|probably woody",
    "^h$|herb"
)

# search for the strings defined before
searchResults <- sapply(searchNames, grepl, oriVals, perl = TRUE)

# name columns of searchResults matrix like woodiness categories
#"non-woody", "woody base", "epiphyte", "hemi-epiphyte", "liana", "graminoid", "woody", "herb"
colnames(searchResults) <- c(1, 2, 1, 1, 1, 1, 3, 1)

# show the number of matches to each category
colSums(searchResults)

# show the original entries for which no match was retrieved
oriVals[rowSums(searchResults) < 1]

# show the number of entries that weren't matched to any category
sum(rowSums(searchResults) < 1)

# show the number of entries that were matched to more that one category
sum(rowSums(searchResults) > 1)

searchResults[rowSums(searchResults) > 1, ] <- FALSE

# use the searchResults matrix to create new value strings by concatenating all data found
newVals <- sapply(seq_len(nrow(searchResults)), function(x) {
    paste(colnames(searchResults)[searchResults[x, ]], collapse = ",")
})
newVals[newVals == ""] <- NA

#Join into TRYSubset
TRYSubset[TraitName == "Plant woodiness", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_woodiness <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_woodiness

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_woodiness, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Leaf type

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Leaf type"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

# define the search strings
searchNames <- c(
    "aphyll|leafless|without leaves",
    "broad\\-?lea(f|v)|^b$|^broad( ?leaved)?$",
    "fine|lanceolate|linear|narrow\\-?lea(f|v)|^n$",
    "needle",
    "scale|micro",
    "conifer|^C$",
    "cylindric",
    "emersed|floating|submersed",
    "photosynthetic stem"
)

# search for the strings defined before
searchResults <- sapply(searchNames, grepl, oriVals)

# name columns of searchResults matrix like corrected searchNames
# "aphyllous", "elliptic", "oblong", "needle", "scale", "coniferous", "cylindric", "aquatic", "photosynthetic stem"
colnames(searchResults) <- c(1, 1, 0, 0, 0, 0, 0, 1, 1)

# show the number of matches to each category
colSums(searchResults)

# show the original entries for which no match was retrieved
oriVals[rowSums(searchResults) < 1]

# show the number of entries that weren't matched to any category
sum(rowSums(searchResults) < 1)

# show the number of entries that were matched to more that one category
sum(rowSums(searchResults) > 1)

searchResults[rowSums(searchResults) > 1, ] <- FALSE

# use the searchResults matrix to create new value strings by concatenating all data found
newVals <- sapply(seq_len(nrow(searchResults)), function(x) {
    paste(colnames(searchResults)[searchResults[x, ]], collapse = ",")
})
newVals[newVals == ""] <- NA

#Join into TRYSubset
TRYSubset[TraitName == "Leaf type", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_leaf <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_leaf

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_leaf, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Lifespan

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Plant lifespan (longevity)"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

searchNames <- c("perennial", "annual", "biennial")
for (i in seq_along(searchNames)) {
    oriVals[grepl(searchNames[i], TRYSubset$OriglName, ignore.case = TRUE) & oriVals == "yes"] <- searchNames[i]
}

# remove ranges and larger signs
oriValsN <- gsub("years?", "", oriVals)
oriValsN <- gsub(".*-", "", oriValsN)
oriValsN <- gsub(">", "", oriValsN)
oriValsN <- gsub("<", "", oriValsN)
oriValsN <- gsub(".*,", "", oriValsN)

# convert strings into numeric values
oriValsN <- as.numeric(oriValsN)

# add information to categorical data
oriVals[oriValsN <= 1] <- "annual" # <= one year
oriVals[oriValsN > 1 & oriValsN <= 2] <- "biennial" #>= one year and <= two years
oriVals[oriValsN > 2 & oriValsN <= 5] <- "perennial" # two years to five years
oriVals[oriValsN > 5] <- "perennial" # five or more years
rm(oriValsN)

# create a vector containing the search strings to look for
searchNames <- c(
    "(^| |\\W)per(e)?(nnial)?|long|woody|tree|shrub|hemicryptophyte|poly-?annual|pluri-?ennial",
    "(^| )an(n)?(ual)?|short|ephemeral",
    "b(i)?-?(asa|e)nnial|bia?s?a?nnual"
)

# search for the strings defined before
searchResults <- sapply(searchNames, grepl, oriVals)

# name columns of searchResults matrix like corrected searchNames
colnames(searchResults) <- c(3, 1, 2)

# show the number of matches to each category
colSums(searchResults)

# show the original entries for which no match was retrieved
oriVals[rowSums(searchResults) < 1]

# show the number of entries that were not matched to any category
sum(rowSums(searchResults) < 1)

# show the number of entries that were matched to more than one category
sum(rowSums(searchResults) > 1)

# use the searchResults matrix to create new value strings by concatenating all data found
newVals <- sapply(seq_len(nrow(searchResults)), function(x) {
    paste(colnames(searchResults)[searchResults[x, ]], collapse = ",")
})
newVals[newVals == ""] <- NA

# Transformar newVals
newVals <- ifelse(grepl("3", newVals), 3, newVals)
newVals <- ifelse(grepl("2", newVals), 2, newVals)

#Join into TRYSubset
TRYSubset[TraitName == "Plant lifespan (longevity)", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_life <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_life

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_life, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Root type

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Root type, root architecture"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

# remove root complexity entries because hard to compare
oriVals[TRYSubset$DatasetID == 115] <- NA

oriVals[TRYSubset$DatasetID == 445 & oriVals == "no"] <- "no cluster roots"
oriVals[TRYSubset$DatasetID == 445 & oriVals == "yes"] <- "cluster roots"

oriVals[TRYSubset$DatasetID == 243 & oriVals == "no"] <- "no taproots"
oriVals[TRYSubset$DatasetID == 243 & oriVals == "yes"] <- "taproots"

oriVals[TRYSubset$DatasetID == 319 & oriVals == "no"] <- "no cluster roots"

# remove purely numeric values and others that have no lowercase character included
oriVals[!grepl("[[:lower:]]", oriVals)] <- NA

# create a vector containing the search strings to look for
searchNames <- c(
    "(^tap ?roots?|root tap|^tap$)",
    "heartroot|intermediate|var\\.root",
    "fibrous",
    "shallow?|lateral",
    "adventitious",
    "^cluster",
    "contractile",
    "no taproots",
    "no cluster roots",
    "absent|no root"
)

# search for the strings defined before
searchResults <- sapply(searchNames, grepl, oriVals)

# name columns of searchResults matrix like corrected searchNames
# "taproot", "heart-root", "fibrous roots", "shallow roots", "adventitious roots", "cluster roots", "contractile roots", "no taproot", "no cluster roots", "no root"
colnames(searchResults) <- c(1, 1, 0, 0, 0, 0, 0, 0, 1, 0)

# show the number of matches to each category
colSums(searchResults)

# show the original entries for which no match was retrieved
oriVals[rowSums(searchResults) < 1]

# show the number of entries that were not matched to any category
sum(rowSums(searchResults) < 1)

# show the number of entries that were matched to more than one category
sum(rowSums(searchResults) > 1)

# use the searchResults matrix to create new value strings by concatenating all data found
newVals <- sapply(seq_len(nrow(searchResults)), function(x) {
    paste(colnames(searchResults)[searchResults[x, ]], collapse = ",")
})
newVals[newVals == ""] <- NA

#Join into TRYSubset
TRYSubset[TraitName == "Root type, root architecture", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_root <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_root

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_root, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Plant life form

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Plant life form (Raunkiaer life form)"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

# repair coded entries
datNames <- names(table(TRYSubset[oriVals == TRUE]$DataName))
for (i in seq_along(datNames)) {
    print(datNames[i])
    print(TRYSubset[DataName == datNames[i] & oriVals == TRUE][1:2])
    print("--------------------------------------")
}
for (i in seq_along(datNames)) {
    oriVals[TRYSubset$DataName == datNames[i] & oriVals == TRUE] <- sub(".*: ", "", datNames[i])
    oriVals[TRYSubset$DataName == datNames[i] & oriVals == "yes"] <- sub(".*: ", "", datNames[i])
}

# remove purely numeric values
oriVals[!grepl("[[:lower:]]", oriVals)] <- NA # remove numeric values

# repair hemi-epiphyte matching problem
oriVals <- gsub("emi-epiph", "emiepiph", oriVals)

newVals <- matrix(NA, length(oriVals), 6)

# get strings with "phyte" ending
unique(tolower(unlist(regmatches(oriVals, gregexpr("\\w+phyte", oriVals)))))

# create a vector containing the search strings to look for
searchNames <- c(
    "chama?ephyte|(^|\\W)cha?(\\W|$)",
    "geophyte|(^|\\W)g(\\W|$)",
    "hemicryptophyte",
    "(?<!pseudo)phanerophyte?|shrub|tree|(^|\\W)ph(\\W|$)",
    "th?erophyte|annual|(^|\\W)th?(\\W|$)"
)

# search for the strings defined before
searchResults <- sapply(searchNames, grepl, oriVals, perl = TRUE)

# name columns of searchResults matrix like corrected searchNames
#   "chamaephyte", "geophyte", "hemicryptophyte", "phanerophyte", "therophyte", "hydrophyte", "cryptophyte", "helophyte"
colnames(searchResults) <- c(4, 2, 3, 5, 1)

# show the number of matches to each category
colSums(searchResults)

# show the original entries for which no match was retrieved
oriVals[rowSums(searchResults) < 1]

# show the number of entries that were not matched to any category
sum(rowSums(searchResults) < 1)

# show the number of entries that were matched to more than one category
sum(rowSums(searchResults) > 1)

# check which entries were classified into > 1 groups
table(oriVals[rowSums(searchResults) > 1])

# remove entries classified into > 1 groups
searchResults[rowSums(searchResults) > 1, ] <- FALSE

# use the searchResults matrix to create new value strings by concatenating all data found
newVals <- sapply(seq_len(nrow(searchResults)), function(x) {
    paste(colnames(searchResults)[searchResults[x, ]], collapse = ",")
})
newVals[newVals == ""] <- NA

#Join into TRYSubset
TRYSubset[TraitName == "Plant life form (Raunkiaer life form)", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_raunkier <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_raunkier

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_raunkier, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Plant growth rate

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Plant growth rate"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

oriVals[!grepl("[[:lower:]]", oriVals)] <- NA

# create a vector containing the search strings to look for
searchNames <- c("low|slow", "medium|moderate", "high|rapid")

# search for the strings defined before
searchResults <- sapply(searchNames, grepl, oriVals)

# name columns of searchResults matrix like corrected searchNames
colnames(searchResults) <- c(1, 2, 3)

# show the number of matches to each category
colSums(searchResults)

# show the original entries for which no match was retrieved
sort(table(oriVals[rowSums(searchResults) < 1]))

# show the number of entries that weren't matched to any category
sum(rowSums(searchResults) < 1)

# show the number of entries that were matched to more that one category
sum(rowSums(searchResults) > 1)

# use the searchResults matrix to create new value strings by concatenating all data found
newVals <- sapply(seq_len(nrow(searchResults)), function(x) {
    paste(colnames(searchResults)[searchResults[x, ]], collapse = ",")
})
newVals[newVals == ""] <- NA

#-------------------------------------------------------------------------------

# Select data of interest
TRYSubset2 <- TRYdata_removed[TraitName == "Plant growth rate"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset2$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

# remove non-numeric values
oriVals[!grepl("\\d", oriVals)] <- NA
oriVals <- as.numeric(oriVals)

oriVals_class <- cut(oriVals, breaks = c(-Inf, 0.3, 1.1, Inf),  # Definir los cortes
                    labels = c(1, 2, 3))  # Etiquetas

#Join into TRYSubset
newVals = coalesce(as.numeric(newVals), as.numeric(oriVals_class))
TRYSubset[TraitName == "Plant growth rate", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_growth <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_growth

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_growth, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Plant flowering time

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Plant reproductive phenology timing (flowering time)"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

# transfer data with OriglName Flowering period: <English month> into
# start of flowering (month) and end of flowering (month)
# enter this data into first two fields per species, set other fields to NA
# convert months into numeric
for (i in seq_along(month.name)) TRYSubset[DatasetID == 516, OriglName := sub(paste0("Flowering period: ", month.name[i]), sprintf("%02d", i), OriglName)]
# sort data to make sure all observations will come in correct order
TRYSubset[, oriOrder := seq_len(nrow(TRYSubset))]
setorder(TRYSubset, DatasetID, SpeciesName, OriglName)
# get first and last consecutive observation month per species
vals <- TRYSubset[DatasetID == 516]$OrigValueStr
vals <- matrix(vals, 12, length(vals) / 12)
frange <- apply(vals, 2, function(x) which(x == TRUE))
fbegin <- vapply(frange, function(x) x[1], 1)
frange <- sapply(frange, diff)
frange <- sapply(frange, paste, collapse = "")
frange <- vapply(frange, nchar, 0)
fend <- sapply(seq_along(fbegin), function(x) fbegin[x] + frange[x])
# replace original data and OriglName
TRYSubset[DatasetID == 516, OrigValueStr := NA]
TRYSubset[DatasetID == 516 & OriglName == "01", OrigValueStr := fbegin]
TRYSubset[DatasetID == 516 & OriglName == "02", OrigValueStr := fend]

# duplicate values with ranges to return start and end separately
TRYTemp <- TRYSubset[OriglName %in% c("Flowering Period", "Flowering time")]
TRYTemp[, OriglName := "temp"]
TRYTemp[, oriOrder := max(TRYSubset$oriOrder) + seq_len(nrow(TRYTemp))]
TRYSubset <- rbind(TRYSubset, TRYTemp)

# create updated oriVals
oriVals <- TRYSubset$OrigValueStr # oriVals == original values
oriVals <- tolower(oriVals)

# new trait names
traitNames <- c(
    "Flowering duration (months)", "Flowering duration (months)",
    "Flowering start (month)", "Flowering start (month)",
    "Flowering end (month)", "Flowering end (month)",
    "Flowering period (season)", "Flowering peak (month)", "Flowering peak (month)",
    "Fruiting start (month)", "Fruiting start (season)", "Fruiting end (season)",
    "Fruiting peak (season)", "Germination (season)"
)

# trait categories based on OriglName
# some are joined later on in case there are different units (days -> months)
origlNames <- c(
    # Flowering duration (months)
    "^(# of Flowering Months|# of Months|FlowDur|flowering duration \\(months\\)|num\\.flring\\.mnths|Reproductive phenology \\(b\\))$",
    # Flowering duration (days)
    "^(Duration of flowering period  \\(visible anthers\\)|Flowering Lasting \\(days\\))$",
    # Flowering start (day)
    "^(Beginning of flowering period  \\(visible anthers\\)|Fdate|onset of flowering|First Flowering Date|Flowering-beginning \\(DOY\\)|FloweringPhenologyStart|FlowStart)$",
    # Flowering start (month)
    "^(01|bloom\\.start|Flower|Flowering|flowering \\(start\\)|Flowering Period|Flowering Period begin|Flowering time|Flowering time: 1\\. earliest month|Onset of Flowering|onset of flowering period|Reproductive phenology \\(a\\))$",
    # Flowering end (day)
    "^(Flowering-End \\(DOY\\))$",
    # Flowering end (month)
    "^(02|end of flowering|FloweringPhenologyEnd|flowering \\(end\\)|Flowering Period end|Flowering time: 2\\. latest month|temp)$",
    # Flowering period (season)
    "^(Bloom Period|Flowering season|FT, flowering time|JAHRESZEIT_E|start of flowering)$",
    # Flowering peak (day)
    "^(FlowerDate|Flowering Period \\(Mediterranean Europe\\)|FlrDate)$",
    # Flowering peak (month)
    "^(Flowering time: 3\\. peak month|Phenology : reproductive)$",
    # Fruiting start (day)
    "^(Smdate)$",
    # Fruiting start (season)
    "^(Fruit/Seed Period Begin)$",
    # Fruiting end (season)
    "^(Fruit/Seed Period End)$",
    # Fruiting peak (season)
    "SE: main seed release period|Time of seed dispersal \\(season\\)",
    # Germination (season)
    "GE, main germination period|Time of germination \\(season\\)"
)

# pre-process data
newVals <- rep(NA, length(oriVals))
for (i in seq_along(origlNames)) {
    # only work on data with the respective OriglName
    tempVals <- oriVals
    tempVals[!grepl(origlNames[i], TRYSubset$OriglName)] <- NA
    # print(sort(table(tempVals)))
    # print(table(TRYSubset[!is.na(tempVals)]$OrigUnitStr))
    if (i == 1) {
        tempVals[tempVals %in% c("/", "na")] <- ""
    } else if (i == 2) {
        tempVals <- as.numeric(format(as.Date(as.numeric(tempVals)), "%m"))
    } else if (i == 3) {
        # convert day to month
        tempVals[!grepl("month", TRYSubset$OrigUnitStr)] <- as.numeric(format(as.Date(as.numeric(
            tempVals[!grepl("month", TRYSubset$OrigUnitStr)]
        )), "%m"))
    } else if (i == 4) {
        # convert names to numbers, ranges to starts
        tempVals[tempVals %in% c("0", "49")] <- ""
        tempVals <- sub("mid-", "", tempVals)
        tempVals <- sub("(\\-|to ).*", "", tempVals)
        tempVals <- sub("^na$", "", tempVals)
        tempVals <- sub("`", "", tempVals)
        temp <- match(tempVals, tolower(month.abb))
        tempVals[!is.na(temp)] <- temp[!is.na(temp)]
        temp <- match(tempVals, tolower(month.name))
        tempVals[!is.na(temp)] <- temp[!is.na(temp)]
        tempVals <- gsub("\\D", "", tempVals)
    } else if (i == 5) {
        # convert day to month
        tempVals[!grepl("month", TRYSubset$OrigUnitStr)] <- as.numeric(format(as.Date(as.numeric(
            tempVals[!grepl("month", TRYSubset$OrigUnitStr)]
        )), "%m"))
    } else if (i == 6) {
        # convert names to numbers, ranges to ends
        tempVals <- sub(".*(\\-|to)", "", tempVals)
        tempVals <- sub("`", "", tempVals)
        temp <- match(tempVals, tolower(month.abb))
        tempVals[!is.na(temp)] <- temp[!is.na(temp)]
        temp <- match(tempVals, tolower(month.name))
        tempVals[!is.na(temp)] <- temp[!is.na(temp)]
        tempVals <- sub("^\\s+", "", tempVals)
    } else if (i == 7) {
        tempVals[tempVals == "not available"] <- ""
        tempVals[tempVals == "jm"] <- "winter"
        tempVals[tempVals == "am"] <- "spring"
        tempVals[tempVals == "ja"] <- "summer"
        tempVals <- sub("any time with rain", "after rain", tempVals)
        tempVals <- sub("all year", "spring, summer, autumn, winter", tempVals)
        tempVals <- sub("autumn-spring", "spring, autumn", tempVals)
        tempVals <- gsub("((start|end) of )?((early|late) )?|mid\\s*|pre-", "", tempVals)
        tempVals <- gsub("fall", "autumn", tempVals)
        tempVals <- gsub("indeterminate|irregular", "variable", tempVals)
        tempVals <- gsub("midsummer", "summer", tempVals)
        tempVals <- gsub("spring-autumn", "spring, summer, autumn", tempVals)
        tempVals <- gsub("-", ", ", tempVals)
    } else if (i == 8) {
        tempVals[tempVals %in% c("", "/")] <- ""
        tempVals[grepl("/", tempVals)] <-
            as.numeric(format(as.Date(tempVals[grepl("/", tempVals)], format = "%d/%m/%Y"), "%j"))
        # convert day to month
        tempVals <- as.numeric(format(as.Date(as.numeric(tempVals)), "%m"))
    } else if (i == 10) {
        # convert day to month
        tempVals <- as.numeric(format(as.Date(as.numeric(tempVals)), "%m"))
    } else if (i %in% c(11:14)) {
        tempVals[tempVals == "not observed"] <- ""
        tempVals[tempVals == "wind"] <- ""
        tempVals <- sub("any time with rain", "after rain", tempVals)
        tempVals <- gsub("fall", "autumn", tempVals)
        tempVals <- sub("year round", "spring, summer, autumn, winter", tempVals)
        tempVals <- sub("no main period", "variable", tempVals)
    }
    newVals[!is.na(tempVals)] <- tempVals[!is.na(tempVals)]
}
newVals[newVals == ""] <- NA

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(newVals)
valueOverview[order(valueOverview)]

# recode text items
newVals <- recode(newVals, "summer" = "7", "spring" = "4", "autumn" = "10", "winter" = "1", "spring, winter" = "", "spring, summer" = "", "autumn, winter" = "")

newVals <- as.numeric(newVals)

#Join into TRYSubset
TRYSubset[TraitName == "Plant reproductive phenology timing (flowering time)", StdValue := newVals]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(StdValue))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, StdValue, UnitName, showOverview = FALSE)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(StdValue), values_fn = list(StdValue = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_flowering <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_flowering

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_flowering, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Litter decomposition rate

# Select data of interest
TRYSubset <- TRYdata_removed[TraitName == "Litter decomposition rate" & OriglName == "k"]

# extract original data strings and change to lowercase to ease classification
oriVals <- TRYSubset$OrigValueStr
oriVals <- tolower(oriVals)

# get an overview over the data by summarizing values and showing them in alphabetical order
valueOverview <- table(oriVals)
valueOverview[order(valueOverview)]

#-------------------------------------------------------------------------------

#Select trait records that have standardized numeric values
num_traits <- rtry_select_row(TRYSubset, complete.cases(TraitID) & complete.cases(OrigValueStr))

#Select the relevant columns for transformation, while suppress the data overview display
num_traits <- rtry_select_col(num_traits, ObservationID, AccSpeciesID, AccSpeciesName, TraitID, TraitName, OrigValueStr, UnitName, showOverview = FALSE)

num_traits$OrigValueStr <- as.numeric(num_traits$OrigValueStr)

#Perform wide table transformation of TraitID, TraitName, and UnitName based on
#ObservationID, AccSpeciesID, and AccSpeciesName with cell values from StdValue
#if several records with StdValue were provided for one trait with the same
#ObservationID, AccSpeciesID, and AccSpeciesName, calculate their mean
num_traits_wider <- rtry_trans_wider(num_traits, names_from = c(TraitID, TraitName, UnitName), values_from = c(OrigValueStr), values_fn = list(OrigValueStr = mean))

#Get trait names
traits<-names(num_traits_wider)

#Get mean values per trait and species
result_decomp <- num_traits_wider %>% 
  group_by(AccSpeciesName) %>% 
  summarise(across(all_of(traits[4:length(traits)]), ~mean(.x, na.rm = TRUE, .names = "{col}")))

result_decomp

#Join trait data to traitxspecies matrix
Trait_matrix_incomplete <- Trait_matrix_incomplete %>%
  full_join(result_decomp, by = "AccSpeciesName")  # "ID" debe estar en ambas tablas

Trait_matrix_incomplete

Data export

Guardamos la matriz de traits x species final.

#Export data
#rtry_export(Trait_matrix_incomplete, file.path('D:/DATOS IPE 2019/Proyecto DIGIPIR/Análisis/TRY', "TRYdata_processed_final.csv"))

Trait matrix processing

---

En este apartado procesamos la matriz de traits x species para rellenar huecos.

Phylogenetic tree

#Create our species list
sp.list <- data.frame(
  species = rownames(Trait_matrix_final),
  genus = NA,
  family = NA,
  species.relative = NA,
  genus.relative = NA
)

#Import plants megatree
megatree <- read.tree('https://raw.githubusercontent.com/megatrees/plant_20221117/main/plant_megatree.tre')

#Import genus relationships
gen.list <- read.csv('https://raw.githubusercontent.com/megatrees/plant_20221117/main/plant_genus_list.csv', sep=",")

#Build phylogenetic tree
result <- phylo.maker(sp.list, megatree, gen.list, nodes.type = 1, scenario = 3)

## [1] "Note: 3 species fail to be binded to the tree."
## [1] "Boraginaceae_sp."        "CAPSELLA_BURSA-PASTORIS"
## [3] "Poaceae_sp"

tree<-result$phylo

#Plot tree
ggtree(tree, layout="circular") +
  geom_tiplab(size=3)

Trait data imputation

#Impute trait values from phylogenetic tree
trait_complete <- impute(Trait_matrix_final, tree, nEigen = 9, addingSpecies = TRUE)

##  adding species 1 / 3 to phylogeny

##  adding species 2 / 3 to phylogeny

##  adding species 3 / 3 to phylogeny

## 
## 0 species added to phylogeny

trait_complete$imputed

#Export data
#rtry_export(trait_complete$imputed, file.path('D:/DATOS IPE 2019/Proyecto DIGIPIR/Análisis/TRY', "TRYdata_imputed_final.csv"))

Trait matrix visualization

---

Complete trait matrix

#Trait Correlation matrix
ggcorr(scale(Trait_matrix_complete_final[c(3:length(Trait_matrix_complete_final))], center=TRUE, scale=TRUE),
  method=c("pairwise.complete.obs", "spearman"),
  high="#3B9AB2", low="#F21A00",
  hjust = 0.75, size=3.5, label=TRUE, label_size=3, label_alpha=TRUE, label_round=1)

#Run PCA
pca<-prcomp(Trait_matrix_complete_final[c(3:length(Trait_matrix_complete_final))], center=TRUE, scale=TRUE)

#PCA plot
autoplot(pca, data = Trait_matrix_complete_final, colour = 'Group', shape = 'Site', size = 3,
         loadings = TRUE, loadings.colour = 'grey25', loadings.label = TRUE, loadings.label.size = 3, 
         loadings.label.colour = 'grey25') +
          theme_bw()

#PCA plot
autoplot(pca, data = Trait_matrix_complete_final, colour = 'grey50', shape = FALSE, label.size = 3,
         loadings = TRUE, loadings.colour = 'blue4', loadings.label = TRUE, loadings.label.size = 3, 
         loadings.label.colour = 'blue4') +
          theme_bw()

#PCA results visual check
fviz_eig(pca, col.var="blue4", addlabels=TRUE)

#Contributions of variables to PC1
fviz_contrib(pca, choice = "var", axes = 1)

#Contributions of variables to PC2
fviz_contrib(pca, choice = "var", axes = 2)

Reduced trait matrix

Tenemos demasiados traits. Muchos de ellos parecen redundantes. Vamos a quedarnos con una selección más pequeña. Los traits que finalmente han sido elegidos son habituales en la bibliografía, se relacionan muy bien con el espectro económico de la hoja (son especies resource-adquisitive or conservative) y están bien representados. Es decir, hemos conseguido información de esos traits para el 70-80% de las especies (aunque de los nutrientes en la hoja -CNP- tenemos información para alrededor del 50% de las especies).

#Trait Correlation matrix
ggcorr(scale(Trait_matrix_complete_final_reduced[c(3:length(Trait_matrix_complete_final_reduced))], center=TRUE, scale=TRUE),
  method=c("pairwise.complete.obs", "spearman"),
  high="#3B9AB2", low="#F21A00",
  hjust = 0.75, size=3.5, label=TRUE, label_size=3, label_alpha=TRUE, label_round=1)

#Run PCA
pca<-prcomp(Trait_matrix_complete_final_reduced[c(3:length(Trait_matrix_complete_final_reduced))], center=TRUE, scale=TRUE)

#PCA plot
autoplot(pca, data = Trait_matrix_complete_final_reduced, colour = 'Group', shape = 'Site', size = 3,
         loadings = TRUE, loadings.colour = 'grey25', loadings.label = TRUE, loadings.label.size = 3, 
         loadings.label.colour = 'grey25') +
          theme_bw()

#PCA plot
autoplot(pca, data = Trait_matrix_complete_final_reduced, colour = 'grey50', shape = FALSE, label.size = 3,
         loadings = TRUE, loadings.colour = 'blue4', loadings.label = TRUE, loadings.label.size = 3, 
         loadings.label.colour = 'blue4') +
          theme_bw()

#PCA results visual check
fviz_eig(pca, col.var="blue4", addlabels=TRUE)

#Contributions of variables to PC1
fviz_contrib(pca, choice = "var", axes = 1)

#Contributions of variables to PC2
fviz_contrib(pca, choice = "var", axes = 2)

Functional diversity

---

Community weighted means (CWMs)

---

En este apartado se calcula el CWM de los valores de los traits para cada plot de 1x1 m2.

Calculation

#Define some parameters
trait_matrix <- Trait_matrix_complete_final_reduced[c(3:length(Trait_matrix_complete_final_reduced))]
sp_matrix <- as.matrix(Plot_sp_matrix_s) #NO PANICULATA!!

#Compute CWMs
CWMs <- functcomp(trait_matrix, sp_matrix, CWM.type = "all")

CWMs

Dimensionality reduction

#Trait Correlation matrix
CWMs_graz <- CWMs %>%
  mutate(Grazing_pressure = Plot_variables_s$Grazing_pressure) %>%
  mutate(Location = Plot_variables_s$Location) %>%
  mutate(Plot_rev = Plot_variables_s$Plot_rev)
ggcorr(scale(CWMs_graz[c(1:(length(CWMs_graz)-2))], center=TRUE, scale=TRUE),
  method=c("pairwise.complete.obs", "spearman"),
  high="#3B9AB2", low="#F21A00",
  hjust = 0.75, size=3.5, label=TRUE, label_size=3, label_alpha=TRUE, label_round=1)

#Run PCA
pca<-prcomp(CWMs, center=TRUE, scale=TRUE)

#PCA plot
autoplot(pca, data = CWMs_graz, size = 3, colour = 'Plot_rev', shape = 'Location',
         loadings = TRUE, loadings.colour = 'grey25', loadings.label = TRUE, loadings.label.size = 3, 
         loadings.label.colour = 'grey25') +
          theme_bw()

#PCA plot
autoplot(pca, data = CWMs_graz, colour = 'grey50', shape = FALSE, label.size = 3,
         loadings = TRUE, loadings.colour = 'blue4', loadings.label = TRUE, loadings.label.size = 3, 
         loadings.label.colour = 'blue4') +
          theme_bw()

#PCA results visual check
fviz_eig(pca, col.var="blue4", addlabels=TRUE)

#Contributions of variables to PC1
fviz_contrib(pca, choice = "var", axes = 1)

#Contributions of variables to PC2
fviz_contrib(pca, choice = "var", axes = 2)

Scores extraction

Podemos resumir los CWMs de cada plot de 1x1 m2 con los valores del primer componente. Valores positivos estarían señalando subcomunidades compuestas por especies resource-adquisitive, caracterizadas por ser de vida más corta, tener hojas más anchas que gruesas, ser algo más pequeñas, menos leñosas, baja relación CN y florecer antes. Las leguminosas tienden a ser de este tipo de especies. Por otra parte, valores negativos estarían señalando subcomunidades dominadas por especies resource-conservative, caracterizadas por ser de vida algo más larta, tener hojas gruesas, ser algo más altas y leñosas, alta relación CN y florecer más tarde.

#Get PC1 values
pca_scores <- data.frame(ID_subplot = Plot_variables_s$ID_subplot, CWMs_PC1 = pca$x[, 1])

#Add scores to the variables data-frame
Plot_variables_s <- Plot_variables_s %>%
  left_join(pca_scores, by = "ID_subplot")

Functional diversity indices

---

En este apartado se calculan diferentes índices de diversidad funcional.

Calculation

#Create a dataframe with trait names and type
traits_t <- data.frame(
  trait_name = colnames(Trait_matrix_complete_final_reduced[c(3:length(Trait_matrix_complete_final_reduced))]),  # trait names
  trait_type = "Q"  # all traits are quantitative
)

#Compute trait-based distances between species
sp_dist <- mFD::funct.dist(
  sp_tr         = Trait_matrix_complete_final_reduced[c(3:length(Trait_matrix_complete_final_reduced))],
  tr_cat        = traits_t,
  metric        = "euclidean",
  scale_euclid  = "scale_center",
  ordinal_var   = "classic",
  weight_type   = "equal",
  stop_if_NA    = TRUE)

#Compute multidimensional trait spaces
fspaces_quality <- mFD::quality.fspaces(
  sp_dist             = sp_dist,
  maxdim_pcoa         = 10,
  deviation_weighting = "absolute",
  fdist_scaling       = FALSE,
  fdendro             = "average")

fspaces_quality$quality_fspaces #Keep 4 dimensions

#I would choose 6. However, it must be fewer than species in plots
#Alternatively, those plots could be removed

#Get PCoA coordinates
sp_faxes_coord <- fspaces_quality$"details_fspaces"$"sp_pc_coord"

#Remove a couple of plots in Cazorla. There must be more species than principal components retained
sp_matrix <- sp_matrix[-c(121, 124), ]

#Compute functional diversity indices in multidimensional trait spaces
alpha_fd_indices <- mFD::alpha.fd.multidim(
  sp_faxes_coord   = sp_faxes_coord[ , c("PC1", "PC2", "PC3", "PC4")],
  asb_sp_w         = sp_matrix,
  ind_vect         = c("fdis", "fmpd", "fnnd", "feve", "fric", "fdiv", "fori", 
                       "fspe", "fide"),
  scaling          = TRUE,
  check_input      = TRUE,
  details_returned = TRUE)

FD indices extraction

#Get FD indices values
fd_ind_values <- alpha_fd_indices$"functional_diversity_indices"

#FD indices correlation matrix
ggcorr(scale(fd_ind_values[c(1,6:7,5,2)], center=TRUE, scale=TRUE),
  method=c("pairwise.complete.obs", "spearman"),
  high="#3B9AB2", low="#F21A00",
  hjust = 0.75, size=3.5, label=TRUE, label_size=3, label_alpha=TRUE, label_round=1)

#Get Fdis values
Fdis <- data.frame(ID_subplot = rownames(fd_ind_values), Fdis = fd_ind_values$fdis)

#Add scores to the variables data-frame
Plot_variables_s <- Plot_variables_s %>%
  left_join(Fdis, by = "ID_subplot")

Exploratory analysis II

---

Data exploration

---

Primary production

Le echamos un vistazo a la relación entre las variables de diversidad funcional y las variables de producción primaria.

#Plot parameters
fig_ndvi_CWM <- ggplot(aes(x = CWMs_PC1, y = NDVI), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,1000) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Dominant traits (CWM)", y ="Primary production (NDVI)", color="Grazing pressure") +
  facet_grid(~Location)

fig_ndvi_Fdis <- ggplot(aes(x = Fdis, y = NDVI), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,1000) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Functional diversity (Fdis)", y ="Primary production (NDVI)", color="Grazing pressure") +
  facet_grid(~Location)

fig_biomass_CWM <- ggplot(aes(x = CWMs_PC1, y = Biomass), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Dominant traits (CWM)", y ="Primary production (biomass)", color="Grazing pressure") +
  facet_grid(~Location)

fig_biomass_Fdis <- ggplot(aes(x = Fdis, y = Biomass), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Functional diversity (Fdis)", y ="Primary production (biomass)", color="Grazing pressure") +
  facet_grid(~Location)

fig_cover_CWM <- ggplot(aes(x = CWMs_PC1, y = Plant_cover), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,100) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Dominant traits (CWM)", y ="Primary production (plant cover)", color="Grazing pressure") +
  facet_grid(~Location)

fig_cover_Fdis <- ggplot(aes(x = Fdis, y = Plant_cover), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  ylim(0,100) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Functional diversity (Fdis)", y ="Primary production (plant cover)", color="Grazing pressure") +
  facet_grid(~Location)

#Print
grid.arrange(fig_ndvi_CWM, fig_ndvi_Fdis, nrow = 2)

grid.arrange(fig_biomass_CWM, fig_biomass_Fdis, nrow = 2)

grid.arrange(fig_cover_CWM, fig_cover_Fdis, nrow = 2)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos, ni 99_mdt_Low_1 y 64_mdt_Low_1 de Cazorla.

Diversity

Le echamos un vistazo a la relación entre las variables de diversidad funcional y las variables de diversidad taxonómica.

fig_CWM_richness <- ggplot(aes(x = Richness, y = CWMs_PC1), data = Plot_variables_s) +
  geom_point(size = 4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Species diversity (richness)", y ="Dominant traits (CWM)", color="Grazing pressure") +
  facet_grid(~Location)

fig_Fdis_richness <- ggplot(aes(x = Richness, y = Fdis), data = Plot_variables_s) +
  geom_point(size = 4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Species diversity (richness)", y ="Functional diversity (Fdis)", color="Grazing pressure") +
  facet_grid(~Location)

fig_CWM_shann <- ggplot(aes(x = EShannon, y = CWMs_PC1), data = Plot_variables_s) +
  geom_point(size = 4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Species diversity (hill numbers)", y ="Dominant traits (CWM)", color="Grazing pressure") +
  facet_grid(~Location)

fig_Fdis_shann <- ggplot(aes(x = EShannon, y = Fdis), data = Plot_variables_s) +
  geom_point(size = 4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Species diversity (hill numbers)", y ="Functional diversity (Fdis)", color="Grazing pressure") +
  facet_grid(~Location)

#Print
grid.arrange(fig_CWM_richness, fig_Fdis_richness, nrow = 2)

grid.arrange(fig_CWM_shann, fig_Fdis_shann, nrow = 2)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos, ni 99_mdt_Low_1 y 64_mdt_Low_1 de Cazorla.

Grazing level

En este apartado se explora la relación entre el nivel de pastoreo y diferentes variables de diversidad funcional.

fig_CWM <- ggplot(aes(x = Plot_rev, y = CWMs_PC1, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Dominant traits (CWM)") +
  facet_grid(~Location)

fig_Fdis <- ggplot(aes(x = Plot_rev, y = Fdis, fill = Plot_rev), data = Plot_variables_s) +
  #geom_violin(trim = FALSE, alpha = .25, colour = NA) +
  geom_jitter(width = .15, aes(color = Plot_rev), alpha = .5, shape = 16, size = 3) +
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = .15, color = "grey25") +
  stat_summary(fun.y = mean, geom = "line", size = .5, aes(group = Location), color = "grey25") +
  stat_summary(fun = mean, geom = "point", size = 4, aes(fill = Plot_rev), shape = 22) + 
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  scale_fill_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Functional diversity (Fdis)") +
  facet_grid(~Location)

#Print
grid.arrange(fig_CWM, fig_Fdis, nrow = 2)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos, ni 99_mdt_Low_1 y 64_mdt_Low_1 de Cazorla.

Grazing level (quantitative)

En este apartado se explora la relación entre el nivel de presión ganadera (como variable cuantitativa) y diferentes variables de diversidad funcional.

#Plot parameters
fig_CWM_q <- ggplot(aes(x = Grazing_pressure, y = CWMs_PC1), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Dominant traits (CWM)", color="Grazing pressure") +
  facet_grid(~Location)

fig_Fdis_q <- ggplot(aes(x = Grazing_pressure, y = Fdis), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Grazing pressure", y ="Functional diversity (Fdis)", color="Grazing pressure") +
  facet_grid(~Location)

#Print
grid.arrange(fig_CWM_q, fig_Fdis_q, nrow = 2)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos, ni 99_mdt_Low_1 y 64_mdt_Low_1 de Cazorla.

Soil fertility

Le echamos un vistazo a la relación entre variables de fertilidad del suelo y variables de diversidad funcional.

fig_CWM_cn <- ggplot(aes(x = CN, y = CWMs_PC1), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (CN)", y ="Dominant traits (CWM)", color="Grazing pressure") +
  facet_grid(~Location)

fig_Fdis_cn <- ggplot(aes(x = CN, y = Fdis), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (CN)", y ="Functional diversity (Fdis)", color="Grazing pressure") +
  facet_grid(~Location)

fig_CWM_p <- ggplot(aes(x = NP, y = CWMs_PC1), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (ln(NP))", y ="Dominant traits (CWM)", color="Grazing pressure") +
  facet_grid(~Location)

fig_Fdis_p <- ggplot(aes(x = NP, y = Fdis), data = Plot_variables_s) +
  geom_point(size =4, aes(color = Plot_rev), alpha = .5, shape = 19) +
  geom_smooth(method="loess", color="black", se=F) +
  geom_smooth(method="lm", se=F) +
  scale_color_manual(values=c("#3aba1e","#a1da27","#f0f123","#f5a527","#db2f14")) +
  theme_bw() +
  theme(strip.background=element_blank()) +
  labs(x ="Soil fertility (ln(NP))", y ="Functional diversity (Fdis)", color="Grazing pressure") +
  facet_grid(~Location)

#Print
grid.arrange(fig_CWM_cn, fig_Fdis_cn, nrow = 2)

grid.arrange(fig_CWM_p, fig_Fdis_p, nrow = 2)

Nota: No se incluyen los plots Low, MedLow y MedHigh de Paniculata en pirineos, ni 99_mdt_Low_1 y 64_mdt_Low_1 de Cazorla.

Statistical analysis

---

Pre-models

---

Species diversity

Ajustamos un modelo para la riqueza de especies.

#Create the random effect factor
Plot_variables_s <- Plot_variables_s %>%
  mutate(Random = paste(Transect, Plot, sep = "_"))

#Set the model
m_sr <- lmer(log(Richness) ~ Grazing_pressure + CN + (1|Random),
             control=lmerControl(optCtrl=list(maxfun=50000)),
             na.action = na.exclude,
             data = Plot_variables_s)

#Model results
summary(m_sr)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Richness) ~ Grazing_pressure + CN + (1 | Random)
##    Data: Plot_variables_s
## Control: lmerControl(optCtrl = list(maxfun = 50000))
## 
## REML criterion at convergence: 99.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.43548 -0.58965  0.04624  0.48878  1.96999 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Random   (Intercept) 0.17989  0.4241  
##  Residual             0.05322  0.2307  
## Number of obs: 140, groups:  Random, 40
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       1.92109    0.17211 69.31282  11.162  < 2e-16 ***
## Grazing_pressure -0.04938    0.07008 34.75468  -0.705 0.485773    
## CN                0.04112    0.01171 79.31180   3.512 0.000738 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Grzng_
## Grzng_prssr  0.006       
## CN          -0.914 -0.031

anova(m_sr, type="III")

r2(m_sr)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.805
##      Marginal R2: 0.146

#Model effects
plot(allEffects(m_sr), multiline=TRUE, ci.style="bars")

#Model assumptions
check_model(m_sr)

CWMs

Ajustamos un modelo para la CWMs.

#Create the random effect factor
Plot_variables_s <- Plot_variables_s %>%
  mutate(Random = paste(Transect, Plot, sep = "_"))

#Set the model
m_cwm <- lmer(CWMs_PC1 ~ Grazing_pressure + Richness + CN + (1|Random),
             control=lmerControl(optCtrl=list(maxfun=50000)),
             na.action = na.exclude,
             data = Plot_variables_s)

#Model results
summary(m_cwm)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CWMs_PC1 ~ Grazing_pressure + Richness + CN + (1 | Random)
##    Data: Plot_variables_s
## Control: lmerControl(optCtrl = list(maxfun = 50000))
## 
## REML criterion at convergence: 550.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.62789 -0.46717 -0.02076  0.44679  2.95280 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Random   (Intercept) 3.673    1.916   
##  Residual             1.502    1.226   
## Number of obs: 140, groups:  Random, 40
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)       -0.15857    0.84285  62.55591  -0.188   0.8514  
## Grazing_pressure   0.56356    0.32233  37.35949   1.748   0.0886 .
## Richness           0.02598    0.03038 135.88317   0.855   0.3939  
## CN                -0.02060    0.05904  87.70062  -0.349   0.7280  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Grzng_ Rchnss
## Grzng_prssr -0.009              
## Richness    -0.216  0.076       
## CN          -0.795 -0.053 -0.296

anova(m_cwm, type="III")

r2(m_cwm)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.727
##      Marginal R2: 0.059

#Model effects
plot(allEffects(m_cwm), multiline=TRUE, ci.style="bars")

#Model assumptions
check_model(m_cwm)

Nota: parece que no hay efecto de la riqueza de especies en el CWMs.

Functional diversity

Ajustamos un modelo para la Fdis.

#Create the random effect factor
Plot_variables_s <- Plot_variables_s %>%
  mutate(Random = paste(Transect, Plot, sep = "_"))

#Set the model
m_fdis <- lmer(Fdis ~ Grazing_pressure + Richness + CN + (1|Random),
             control=lmerControl(optCtrl=list(maxfun=50000)),
             na.action = na.exclude,
             data = Plot_variables_s)

#Model results
summary(m_fdis)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Fdis ~ Grazing_pressure + Richness + CN + (1 | Random)
##    Data: Plot_variables_s
## Control: lmerControl(optCtrl = list(maxfun = 50000))
## 
## REML criterion at convergence: -377.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.79073 -0.41818  0.03441  0.51509  2.83791 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Random   (Intercept) 0.002968 0.05448 
##  Residual             0.001691 0.04112 
## Number of obs: 138, groups:  Random, 40
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)      1.810e-01  2.542e-02 5.602e+01   7.122 2.17e-09 ***
## Grazing_pressure 7.772e-03  9.361e-03 3.547e+01   0.830    0.412    
## Richness         4.011e-05  9.813e-04 1.330e+02   0.041    0.967    
## CN               9.945e-04  1.813e-03 7.866e+01   0.549    0.585    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Grzng_ Rchnss
## Grzng_prssr -0.006              
## Richness    -0.212  0.084       
## CN          -0.783 -0.061 -0.334

anova(m_fdis, type="III")

r2(m_fdis)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.644
##      Marginal R2: 0.018

#Model effects
plot(allEffects(m_fdis), multiline=TRUE, ci.style="bars")

#Model assumptions
check_model(m_fdis)

Soil fertility

Ajustamos un modelo para la relación CN.

#Create the random effect factor
Plot_variables_s <- Plot_variables_s %>%
  mutate(Random = paste(Transect, Plot, sep = "_"))

#Set the model
m_cn <- lmer(CN ~ Grazing_pressure + (1|Random),
             control=lmerControl(optCtrl=list(maxfun=50000)),
             na.action = na.exclude,
             data = Plot_variables_s)

#Model results
summary(m_cn)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CN ~ Grazing_pressure + (1 | Random)
##    Data: Plot_variables_s
## Control: lmerControl(optCtrl = list(maxfun = 50000))
## 
## REML criterion at convergence: 603.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7002 -0.1752 -0.0351  0.2404  5.8218 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Random   (Intercept) 23.575   4.855   
##  Residual              1.407   1.186   
## Number of obs: 140, groups:  Random, 40
## 
## Fixed effects:
##                  Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)       13.4093     0.7755 37.9190  17.290   <2e-16 ***
## Grazing_pressure   0.1922     0.7763 37.9096   0.248    0.806    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Grzng_prssr -0.055

anova(m_cn, type="III")

r2(m_cn)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.944
##      Marginal R2: 0.001

#Model effects
plot(allEffects(m_cn), multiline=TRUE, ci.style="bars")

#Model assumptions
check_model(m_cn)

Primary production

Ajustamos un modelo para la producción primaria.

#Create the random effect factor
Plot_variables_s <- Plot_variables_s %>%
  mutate(Random = paste(Transect, Plot, sep = "_"))

#Set the model
m_ndvi <- lmer(NDVI ~ Grazing_pressure + Richness + CWMs_PC1 + Fdis + CN + (1|Random),
             control=lmerControl(optCtrl=list(maxfun=50000)),
             na.action = na.exclude,
             data = Plot_variables_s)

#Model results
summary(m_ndvi)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: NDVI ~ Grazing_pressure + Richness + CWMs_PC1 + Fdis + CN + (1 |  
##     Random)
##    Data: Plot_variables_s
## Control: lmerControl(optCtrl = list(maxfun = 50000))
## 
## REML criterion at convergence: 1419.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7944 -0.4709 -0.0701  0.4727  1.9498 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Random   (Intercept) 19590    140.0   
##  Residual             10436    102.2   
## Number of obs: 116, groups:  Random, 34
## 
## Fixed effects:
##                  Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)       408.042     83.791  54.945   4.870 9.85e-06 ***
## Grazing_pressure  -13.120     25.891  16.263  -0.507   0.6191    
## Richness            5.251      2.729 108.648   1.924   0.0569 .  
## CWMs_PC1           14.695      8.692 109.589   1.691   0.0937 .  
## Fdis              331.904    275.182  94.700   1.206   0.2308    
## CN                  8.666      5.062  47.708   1.712   0.0934 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Grzng_ Rchnss CWM_PC Fdis  
## Grzng_prssr -0.042                            
## Richness    -0.163  0.052                     
## CWMs_PC1     0.326 -0.129 -0.076              
## Fdis        -0.573  0.050  0.041 -0.525       
## CN          -0.565 -0.063 -0.364  0.070 -0.129

anova(m_ndvi, type="III")

r2(m_ndvi)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.721
##      Marginal R2: 0.196

#Model effects
plot(allEffects(m_ndvi), multiline=TRUE, ci.style="bars")

#Model assumptions
check_model(m_ndvi)

Pre-models by location

---

Mediante una app interactiva (shiny) ejecutamos los diferentes modelos para cada sitio por separado.

#User interface
ui <- fluidPage(
  titlePanel("LMMs pre-analysis"),
  sidebarLayout(
    sidebarPanel(
      selectInput("location", "Select location:", 
                  choices = unique(Plot_variables_s$Location),
                  selected = "Pirineos"),
      
      selectInput("selected_model", "Select response variable:", 
                  choices = c("Richness", "CWMs_PC1", "Fdis", "CN", "NDVI"),
                  selected = "Richness"),
      
      actionButton("run_model", "Run model")
    ),
    
    mainPanel(
      tabsetPanel(
        tabPanel("Summary", verbatimTextOutput("summary")),
        tabPanel("ANOVA", verbatimTextOutput("anova")),
        tabPanel("R²", verbatimTextOutput("r2")),
        tabPanel("Effects", plotOutput("effects_plot")),
        tabPanel("Assumptions", plotOutput("check_plot"))
      )
    )
  )
)

#Server
server <- function(input, output, session) {
  
  #Create filter data by location
  datos_filtrados <- reactive({
    req(input$location)
    Plot_variables_s %>%
      filter(Location == input$location) %>%
      mutate(Random = paste(Transect, Plot, sep = "_"))
  })
  
  #Define the different models
  modelo <- eventReactive(input$run_model, {
    req(datos_filtrados())
    
    formula <- switch(input$selected_model,
                      "Richness" = log(Richness) ~ Grazing_pressure + CN + (1|Random),
                      "CWMs_PC1" = CWMs_PC1 ~ Grazing_pressure + Richness + CN + (1|Random),
                      "Fdis" = Fdis ~ Grazing_pressure + Richness + CN + (1|Random),
                      "CN" = CN ~ Grazing_pressure + (1|Random),
                      "NDVI" = NDVI ~ Grazing_pressure + Richness + CWMs_PC1 + Fdis + CN + (1|Random))
    
    lmer(formula, 
         control = lmerControl(optCtrl = list(maxfun = 50000)),
         na.action = na.exclude,
         data = datos_filtrados())
  })
  
  #Show results
  output$summary <- renderPrint({ req(modelo()); summary(modelo()) })
  output$anova <- renderPrint({ req(modelo()); anova(modelo(), type = "III") })
  output$r2 <- renderPrint({ req(modelo()); r2(modelo()) })
  output$effects_plot <- renderPlot({ req(modelo()); plot(allEffects(modelo()), multiline = TRUE, ci.style = "bars") })
  output$check_plot <- renderPlot({ req(modelo()); check_model(modelo()) })
}

# Ejecutar la app
shinyApp(ui = ui, server = server)

SEM

---

Sites combined

Hacemos un SEM con todos los datos.

#Create the random effect factor
Plot_variables_s <- Plot_variables_s %>%
  mutate(Random = paste(Transect, Plot, sep = "_"))

#Sem build
sem_all <- psem(
    lmer(log(Richness) ~ Grazing_pressure + CN + (1|Random),
             control=lmerControl(optCtrl=list(maxfun=50000)),
             na.action = na.exclude,
             data = Plot_variables_s),
    lmer(CWMs_PC1 ~ Grazing_pressure + CN + (1|Random),
             control=lmerControl(optCtrl=list(maxfun=50000)),
             na.action = na.exclude,
             data = Plot_variables_s),
    lmer(Fdis ~ Grazing_pressure + CN + (1|Random),
             control=lmerControl(optCtrl=list(maxfun=50000)),
             na.action = na.exclude,
             data = Plot_variables_s),
    lmer(CN ~ Grazing_pressure + (1|Random),
             control=lmerControl(optCtrl=list(maxfun=50000)),
             na.action = na.exclude,
             data = Plot_variables_s),
    lmer(NDVI ~ Grazing_pressure + log(Richness) + CWMs_PC1 + Fdis + CN + (1|Random),
             control=lmerControl(optCtrl=list(maxfun=50000)),
             na.action = na.exclude,
             data = Plot_variables_s),
    data = Plot_variables_s)

#Sem results
summary(sem_all, .progressBar = FALSE)

## 
## Structural Equation Model of sem_all 
## 
## Call:
##   log(Richness) ~ Grazing_pressure + CN
##   CWMs_PC1 ~ Grazing_pressure + CN
##   Fdis ~ Grazing_pressure + CN
##   CN ~ Grazing_pressure
##   NDVI ~ Grazing_pressure + log(Richness) + CWMs_PC1 + Fdis + CN
## 
##     AIC
##  2324.554
## 
## ---
## Tests of directed separation:
## 
##                   Independ.Claim Test.Type       DF Crit.Value P.Value    
##   CWMs_PC1 ~ log(Richness) + ...      coef 133.3687     1.3426  0.1817    
##       Fdis ~ log(Richness) + ...      coef 126.1179     0.5168  0.6062    
##            Fdis ~ CWMs_PC1 + ...      coef 120.7347     7.6113  0.0000 ***
## 
## --
## Global goodness-of-fit:
## 
## Chi-Squared = 31.155 with P-value = 0 and on 3 degrees of freedom
## Fisher's C = 55.887 with P-value = 0 and on 6 degrees of freedom
## 
## ---
## Coefficients:
## 
##        Response        Predictor Estimate Std.Error       DF Crit.Value P.Value
##   log(Richness) Grazing_pressure  -0.0494    0.0701  34.7547    -0.7046  0.4858
##   log(Richness)               CN   0.0411    0.0117  79.3118     3.5117  0.0007
##        CWMs_PC1 Grazing_pressure   0.5427    0.3208  37.1105     1.6914  0.0991
##        CWMs_PC1               CN  -0.0057    0.0563  71.3808    -0.1007  0.9200
##            Fdis Grazing_pressure   0.0077    0.0093  36.7423     0.8342  0.4096
##            Fdis               CN   0.0010    0.0017  63.6929     0.6000  0.5506
##              CN Grazing_pressure   0.1922    0.7763  37.9096     0.2476  0.8058
##            NDVI Grazing_pressure -14.1980   20.7906  12.3515    -0.6829  0.5073
##            NDVI    log(Richness) 139.6031   36.4530  62.2374     3.8297  0.0003
##            NDVI         CWMs_PC1  17.0156    8.5079  99.1059     2.0000  0.0482
##            NDVI             Fdis 291.2328  276.5861 109.0528     1.0530  0.2947
##            NDVI               CN   8.0597    4.5055  30.9518     1.7888  0.0834
##   Std.Estimate    
##        -0.0922    
##         0.3631 ***
##         0.2353    
##        -0.0116    
##         0.1165    
##         0.0727    
##         0.0406    
##        -0.0671    
##         0.3537 ***
##         0.1856   *
##         0.0915    
##         0.1804    
## 
##   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05
## 
## ---
## Individual R-squared:
## 
##   Response method Marginal Conditional
##   Richness   none     0.15        0.81
##   CWMs_PC1   none     0.05        0.73
##       Fdis   none     0.02        0.64
##         CN   none     0.00        0.94
##       NDVI   none     0.37        0.68

#Sem plot
plot(sem_all)

Pirineos

Cazorla